{
"cells": [
{
"cell_type": "code",
"execution_count": 1,
"id": "bab3bf3f-fa2a-4402-8aae-ef8256fedc9d",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"%matplotlib inline\n",
"%config InlineBackend.figure_format = 'svg'\n",
"\n",
"import numpy as np\n",
"from triqs.operators import Operator, n\n",
"from triqs.gf import MeshImFreq, MeshImTime, Gf, make_gf_from_fourier, iOmega_n, inverse, Fourier\n",
"from triqs.plot.mpl_interface import oplot, oploti, oplotr, plt"
]
},
{
"cell_type": "markdown",
"id": "acf8653d-37dd-4636-b45b-b177cef08f6f",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"# One-boson retarded interaction\n",
"\n",
"To showcase the ability to treat retarded interactions in **TRIQS/ctseg**, we consider the minimal discrete model with one fermion (0) coupled to a fermionic bath state (1) as well as a boson\n",
"\n",
"$$\n",
"H = \\epsilon_0 c^\\dagger_0 c_0 + \\epsilon_1 c^\\dagger_1 c_1 \n",
"+ V (c^\\dagger_1 c_0 + c^\\dagger_0 c_1)\n",
"+ g (b + b^\\dagger)c^\\dagger_0 c_0\n",
"+ \\omega_0 b^\\dagger b\n",
"\\, .\n",
"$$\n",
"\n",
"Integrating out the bath fermion and the boson gives the single fermion action \n",
"$$\n",
"\\mathcal{S} = \n",
"\\int_0^\\beta d\\tau H_{0}(\\tau) \n",
"+ \n",
"\\iint_0^\\beta d\\tau d\\tau' c_0^\\dagger(\\tau) \\Delta(\\tau - \\tau') c_0(\\tau)\n",
"+\n",
"\\iint_0^\\beta d\\tau d\\tau' n_0^\\dagger(\\tau) \\mathcal{U}(\\tau - \\tau') n_0(\\tau)\n",
"$$\n",
"with $H_{0} = \\epsilon_0 c^\\dagger_0 c_0$, retarded hybridization \n",
"$\\Delta(i\\nu_n) = \\frac{V^2}{i\\nu_n - \\epsilon_1}$,\n",
"and retarded interaction $\\mathcal{U}(i\\omega_n) = - 2 g^2 \\frac{\\omega_0}{\\omega_0^2 - (i\\omega_n)^2} \\,.$\n",
"\n",
"The fermionic Green's function $G_0(\\tau) = - \\langle \\mathcal{T} c_0(\\tau) c_0\\rangle$ of the action $\\mathcal{S}$ can be computed numerically exactly using **TRIQS/ctseg**, and (since the number of states is finite) also using exact diagonalization of $H$. "
]
},
{
"cell_type": "markdown",
"id": "8d15b5db-3a8a-49ba-a141-15a32c73c6a0",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## Simulating the retarded action $\\mathcal{S}$\n",
"\n",
"To compute the fermionic Green's function $G_0(\\tau)$ with **TRIQS/ctseg** we\n",
"\n",
"**1.** set up the system parameters, $\\epsilon_0$, $\\epsilon_1$, $V$, $g$, and $\\omega_0$,"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "c904ec9d-c533-42ec-afe9-4657d4e9f213",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"eps0, eps1, V, g, omega0 = -0.1, +0.1, 0.25, 1., 1.\n",
"mu = -g**2 / omega0 # For half-filling at eps0 = 0"
]
},
{
"cell_type": "markdown",
"id": "cbb15884-fdf0-4c08-a57a-e6711541ad9d",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"**2.** construct a solver instance,"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "710e79a6-e2b5-4eac-a7d9-1f33a329ebdb",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Starting run with 1 MPI rank(s) at : 2026-03-11 10:21:51.593224\n"
]
}
],
"source": [
"from triqs_ctseg import Solver\n",
"\n",
"S = Solver(beta=10., n_tau_bosonic=1001,\n",
" gf_struct = [[\"0\", 1], [\"1\", 1]]); # CTSEG requires two flavours or more, adding one inactive fermion"
]
},
{
"cell_type": "markdown",
"id": "840d4260-08e9-4966-a74a-a52eeeb1a7a9",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"**3.** set the hybridization function $\\Delta(i\\nu_n) = \\frac{V^2}{i\\nu_n - \\epsilon_1}$,"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "ffad23a2-a375-44f7-8e81-e65f22ffd541",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"Delta_iw = make_gf_from_fourier(S.Delta_tau['0'])\n",
"Delta_iw << V**2 * inverse(iOmega_n - eps1)\n",
"S.Delta_tau['0'] << Fourier(Delta_iw);"
]
},
{
"cell_type": "markdown",
"id": "ab4fa9b7-aaf5-4d54-9ee0-efdafc69a9c6",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"**4.** set the retarded interaction $\\mathcal{U}(i\\omega_n) = - 2g^2 \\frac{\\omega_0}{ \\omega_0^2 - (i\\omega_n)^2}$,"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "29938c11-19b5-4820-9a1b-e0d1b2811df3",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"D0_iw = make_gf_from_fourier(S.D0_tau['0','0'])\n",
"D0_iw << -2 * g**2 * omega0 * inverse(omega0**2 - iOmega_n*iOmega_n);\n",
"S.D0_tau[\"0\", \"0\"] << Fourier(D0_iw);"
]
},
{
"cell_type": "markdown",
"id": "e5ff0fb9-56e7-46b7-b987-40690601fc1c",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"**5.** and finally, run the solver while passing the Hamiltonian $H_{0} = (\\epsilon_0 - \\mu) c^\\dagger_0 c_0$."
]
},
{
"cell_type": "code",
"execution_count": 6,
"id": "0adc4c09-1375-4993-ae2f-6f71dba85240",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\n",
"╔╦╗╦═╗╦╔═╗ ╔═╗ ┌─┐┌┬┐┌─┐┌─┐┌─┐\n",
" ║ ╠╦╝║║═╬╗╚═╗ │ │ └─┐├┤ │ ┬\n",
" ╩ ╩╚═╩╚═╝╚╚═╝ └─┘ ┴ └─┘└─┘└─┘\n",
"\n",
"\n",
"Block: 0 Index: 0 Color: 0\n",
"Block: 1 Index: 0 Color: 1\n",
"\n",
" Interaction matrix: U = \n",
"[[0,0]\n",
" [0,0]] \n",
"\n",
"Orbital energies: mu - eps = [-0.9,-0] \n",
"\n",
"\n",
" Renormalized interaction matrix: U = \n",
"[[0,0]\n",
" [0,0]] \n",
"\n",
"Renormalized orbital energies: mu - eps = [0.100008,0] \n",
"\n",
"Dynamical interactions = true, Jperp interactions = false \n",
"\n",
"[Rank 0] Performing warmup phase...\n",
"[Rank 0] 10:21:51 3.77% done, ETA 00:00:02, cycle 37653 of 1000000, 3.77e+05 cycles/sec\n",
"[Rank 0] 10:21:53 82.59% done, ETA 00:00:00, cycle 825938 of 1000000, 3.89e+05 cycles/sec\n",
"[Rank 0] 10:21:54 100.00% done, ETA 00:00:00, cycle 1000000 of 1000000, 3.90e+05 cycles/sec\n",
"[Rank 0] Performing accumulation phase...\n",
"[Rank 0] 10:21:54 0.30% done, ETA 00:00:33, cycle 30026 of 10000000, 3.00e+05 cycles/sec\n",
"[Rank 0] 10:21:56 6.26% done, ETA 00:00:31, cycle 626337 of 10000000, 2.95e+05 cycles/sec\n",
"[Rank 0] 10:21:58 13.66% done, ETA 00:00:29, cycle 1366104 of 10000000, 2.93e+05 cycles/sec\n",
"[Rank 0] 10:22:02 23.04% done, ETA 00:00:26, cycle 2303792 of 10000000, 2.95e+05 cycles/sec\n",
"[Rank 0] 10:22:05 34.70% done, ETA 00:00:22, cycle 3470255 of 10000000, 2.95e+05 cycles/sec\n",
"[Rank 0] 10:22:10 49.55% done, ETA 00:00:17, cycle 4955324 of 10000000, 2.96e+05 cycles/sec\n",
"[Rank 0] 10:22:17 67.81% done, ETA 00:00:10, cycle 6780908 of 10000000, 2.96e+05 cycles/sec\n",
"[Rank 0] 10:22:24 90.58% done, ETA 00:00:03, cycle 9057730 of 10000000, 2.96e+05 cycles/sec\n",
"[Rank 0] 10:22:27 100.00% done, ETA 00:00:00, cycle 10000000 of 10000000, 2.96e+05 cycles/sec\n",
"[Rank 0] Collect results: Waiting for all MPI processes to finish accumulating...\n",
"Densities:\n",
" 0: [0.700541]\n",
" 1: [0.49906]\n",
"\n",
"[Rank 0] Warmup duration: 2.5641 seconds [00:00:02]\n",
"[Rank 0] Simulation duration: 33.7571 seconds [00:00:33]\n",
"[Rank 0] Number of measures: 10000000\n",
"[Rank 0] Cycles (measures) / second: 2.96e+05\n",
"[Rank 0] Measurement durations (total = 7.3851):\n",
"[Rank 0] Measure : Duration = 6.4573\n",
"[Rank 0] Measure Average Sign: Duration = 0.1683\n",
"[Rank 0] Measure Densities: Duration = 0.1763\n",
"[Rank 0] Measure G(tau)/F(tau): Duration = 0.4049\n",
"[Rank 0] Measure Perturbation order Delta: Duration = 0.1783\n",
"[Rank 0] Move statistics:\n",
"[Rank 0] Move insert: Proposed = 22858922, Accepted = 2849735, Rate = 0.1247\n",
"[Rank 0] Move remove: Proposed = 22852384, Accepted = 2849929, Rate = 0.1247\n",
"[Rank 0] Move double insert: Proposed = 22868112, Accepted = 0, Rate = 0.0000\n",
"[Rank 0] Move double remove: Proposed = 22857446, Accepted = 0, Rate = 0.0000\n",
"[Rank 0] Move move: Proposed = 22852998, Accepted = 1564818, Rate = 0.0685\n",
"[Rank 0] Move split: Proposed = 22858879, Accepted = 4473405, Rate = 0.1957\n",
"[Rank 0] Move regroup: Proposed = 22851259, Accepted = 4473212, Rate = 0.1958\n",
"[Rank 0] Move durations (total = 18.2591):\n",
"[Rank 0] Move insert: Duration = 3.4120\n",
"[Rank 0] Move remove: Duration = 1.9373\n",
"[Rank 0] Move double insert: Duration = 3.4507\n",
"[Rank 0] Move double remove: Duration = 0.9989\n",
"[Rank 0] Move move: Duration = 2.4168\n",
"[Rank 0] Move split: Duration = 3.8561\n",
"[Rank 0] Move regroup: Duration = 2.1873\n",
"\n",
"Total number of measures: 10000000\n",
"Total cycles (measures) / second: 2.96e+05\n",
"Average sign: 1\n",
"Average perturbation order in Delta: 1.523\n"
]
}
],
"source": [
"S.solve(h_loc0=(eps0 - mu) * n(\"0\", 0), h_int=Operator(), # No static interaction\n",
" measure_nn_tau=True, length_cycle=16, n_warmup_cycles=int(1e6), n_cycles=int(1e7))"
]
},
{
"cell_type": "markdown",
"id": "c9297f6b-d837-4b3a-bb20-e88840cc4507",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## Reference solution from exact diagonalization (using **pyED**)\n",
"\n",
"To get the exact solution we setup the full Hamiltonian $H$ with two fermions and one boson, diagonalize $H$, and compute the single particle Green's function $G_0(\\tau)$."
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "0b4a11fe-49c8-4e9d-ad4f-87c7f67c2869",
"metadata": {
"editable": true,
"nbsphinx": "hidden",
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [],
"source": [
"from pyed.SparseExactDiagonalization import SparseExactDiagonalization\n",
"from pyed.SparseMatrixFockStates import SparseMatrixFermiBoseCreationOperators\n",
"\n",
"def get_ed_ref(eps0, eps1, V, g, omega0, mesh_f_tau, mesh_b_tau, Nb_max=10):\n",
"\n",
" ops = SparseMatrixFermiBoseCreationOperators(Nf=2, Nb=1, Nb_max=Nb_max)\n",
" \n",
" c0, c1, b = ops.c_dag[0].getH(), ops.c_dag[1].getH(), ops.b_dag[0].getH()\n",
" n0, n1, nb = c0.getH() * c0, c1.getH() * c1, b.getH() * b\n",
"\n",
" H = eps0 * n0 + eps1 * n1 + V*(c1.getH() * c0 + c0.getH() * c1) + g*(b + b.getH())*n0 + omega0 * nb\n",
"\n",
" ed = SparseExactDiagonalization(H, mesh_f_tau.beta)\n",
"\n",
" tau_f = np.array([float(t) for t in mesh_f_tau])\n",
" g_tau = Gf(mesh=mesh_f_tau, target_shape=[1, 1])\n",
" g_tau.data[:, 0, 0] = ed.get_tau_greens_function_component(tau_f, c0, c0.getH())\n",
"\n",
" tau_b = np.array([float(t) for t in mesh_b_tau])\n",
" chi_tau = Gf(mesh=mesh_b_tau, target_shape=[1, 1])\n",
" chi_tau.data[:, 0, 0] = -ed.get_tau_greens_function_component(tau_b, n0, n0)\n",
" \n",
" return g_tau, chi_tau"
]
},
{
"cell_type": "markdown",
"id": "1e1bec02-660b-43f3-9b63-0a0df9bf7ad6",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"\n",
"\n",
"## Single particle Green's function\n",
"\n",
"Comparing the single particle Green's function computed with **TRIQS/ctseg** and **pyED** we find excellent agreement."
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "850531e2-9236-4536-90da-e7476608ea0d",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"g_tau_ed, chi_tau_ed = get_ed_ref(eps0 - mu, eps1, V, g, omega0, \n",
" S.results.G_tau[\"0\"].mesh, S.results.nn_tau[\"0\", \"0\"].mesh)\n",
"\n",
"oplotr(S.results.G_tau[\"0\"], lw=0.5, alpha=0.75, label=f'TRIQS/ctseg')\n",
"oplotr(g_tau_ed, '--', alpha=0.75, label=f'pyED')\n",
"plt.xlabel(r'$\\tau$'); plt.ylabel(r'$G_0(\\tau)$'); plt.ylim(top=0);"
]
},
{
"cell_type": "markdown",
"id": "cc78f225-5fba-4a55-a233-9763e20ff79e",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"source": [
"## Dynamic density correlation\n",
"\n",
"Also the density correlation function $\\chi(\\tau) = \\langle \\mathcal{T} n_0(\\tau) n_0 \\rangle$ matches the reference solution from exact diagonalization."
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "135c48e8-b585-42dd-b036-3670de78d93e",
"metadata": {
"editable": true,
"slideshow": {
"slide_type": ""
},
"tags": []
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"oplotr(S.results.nn_tau['0', '0'], alpha=0.75, label='TRIQS/ctseg')\n",
"oplotr(chi_tau_ed, '--', alpha=0.75, label='pyED')\n",
"plt.xlabel(r'$\\tau$'); plt.ylabel(r'$\\chi(\\tau)$');"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.12.12"
}
},
"nbformat": 4,
"nbformat_minor": 5
}